module sphere_geometry_mod

	use const_mod

	implicit none

	private
	public calc_distance
  public cartesian_transform
  public rotate2regular
  public regular2rotate

contains

	pure real(8) function calc_distance(lon1, lat1, lon2, lat2) result(res)

		real(8), intent(in) :: lon1
		real(8), intent(in) :: lat1
		real(8), intent(in) :: lon2
		real(8), intent(in) :: lat2

		res = radius * acos(min(1.0d0, max(-1.d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))

	end function calc_distance

	subroutine cartesian_transform(lon, lat, x, y, z)

		real(8), intent(in) :: lon, lat
		real(8), intent(out) :: x, y, z

		real(8) cos_lat

		cos_lat = cos(lat)
		x = radius * cos_lat * cos(lon)
		y = radius * cos_lat * sin(lon)
		z = radius * sin(lat)

	end subroutine cartesian_transform

	subroutine regular2rotate(lon_reg, lat_reg, lonp_rot, latp_rot, lon_rot, lat_rot)

		real(8), intent(in) :: lon_reg, lat_reg
		real(8), intent(in) :: lonp_rot, latp_rot
		real(8), intent(out):: lon_rot, lat_rot

		lon_rot = atan2(cos(lat_reg) * sin(lon_reg - lonp_rot), &
			        cos(lat_reg) * sin(latp_rot) * cos(lon_reg - lonp_rot) - cos(latp_rot) * sin(lat_reg))
		if (lon_rot < 0.0) lon_rot = pi2 + lon_rot

		lat_rot = asin(min(max(sin(lat_reg) * sin(latp_rot) + &
			        cos(lat_reg) * cos(latp_rot) * cos(lon_reg - lonp_rot), -1.0), 1.0))

	end subroutine regular2rotate

	subroutine rotate2regular(lon_rot, lat_rot, lonp_rot, latp_rot, lon_reg, lat_reg)

		real(8), intent(in) :: lon_rot, lat_rot
		real(8), intent(in) :: lonp_rot, latp_rot
		real(8), intent(out):: lon_reg, lat_reg

		lon_reg = lonp_rot + atan2(cos(lat_rot) * sin(lon_rot), &
				                         sin(lat_rot) * cos(latp_rot) + cos(lat_rot) * cos(lon_rot) * sin(latp_rot))
		if (lon_reg > pi2) lon_reg = lon_reg - pi2

		lat_reg = asin(min(max(sin(lat_rot) * sin(latp_rot) -&
				                 cos(lat_rot) * cos(latp_rot) * cos(lon_rot), -1.0), 1.0))

	end subroutine rotate2regular

end module sphere_geometry_mod